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ABSTRACT 

The  dot  diflPusion  method  for  digital  halftoning  has  the 
advantage  of  parallelism  unlike  the  error  diffusion  method. 
However,  image  quality  offered  by  error  diffusion  is  still  re¬ 
garded  as  superior  to  other  known  methods.  In  this  paper  we 
show  how  the  dot  diffusion  method  can  be  improved  by  opti¬ 
mization  of  the  so-called  class  matrix.  By  taking  the  human 
visual  characteristics  into  account  we  show  that  such  opti¬ 
mization  consistently  results  in  images  comparable  to  error 
diffusion,  without  sacrificing  the  parallelism. 

1  INTRODUCTION 

Digital  halftoning  is  the  rendition  of  continuous-tone  pictures 
on  displays  that  are  capable  of  producing  only  two  levels. 
There  are  many  good  methods  for  digital  halftoning:  ordered 
dither  [1],  error  diffusion  [2],  neural-net  based  methods  [8], 
and  more  recently  direct  binary  search  (DBS)  [7].  Ordered 
dithering  is  a  thresholding  of  the  continuous-tone  image  with 
a  spatially  periodic  screen  [1]  .  In  error  diffusion  [2] ,  the  error 
is  'diffused’  to  the  unprocessed  neighbor  points. 

Ordered  dithering  is  a  parallel  method,  requiring  only 
pointwise  comparisons.  But  the  resulting  halftones  suffer 
from  periodic  patterns.  On  the  other  hand  error  diffused 
halftones  do  not  suffer  from  periodicity  and  offer  blue  noise 
characteristic  [3]  which  is  found  to  be  desirable.  The  main 
drawback  is  that  error  diffusion  is  inherently  serial,  i.e.,  to 
get  the  halftoned  value  of  the  last  pixel,  all  of  the  remaining 
points  should  be  processed.  Also  there  occur  worm-like  pat¬ 
terns  in  near  mid-gray  regions  and  resulting  halftones  have 
ghosting  problem  [6].  Mitsa  and  Parker  have  optimized  or¬ 
dered  dither  matrix  [4]  for  large  size  like  256x256  to  get  the 
blue  noise  effect.  This  is  a  compromise  between  parallelism 
and  image  quality. 

The  dot  diffusion  method  for  halftoning,  introduced  by 
Knuth  [6]  is  an  attractive  method  which  attempts  to  retain 
the  good  features  of  error  diffusion  while  offering  substantial 
parallelism.  However,  surprisingly,  not  much  work  has  been 
done  on  optimization  of  the  so-called  class  matrix.  In  this 
work  we  will  show  that  the  class  matrix  (see  below)  can 
further  be  optimized  by  taking  into  account  the  properties 
of  human  visual  system  (HVS).  The  resulting  halftones  will 
then  be  of  the  same  quality  as  for  error  diffusion.  Since  dot 
diffusion  also  offers  increased  pctrallelism,  it  now  appears  to 
be  an  attractive  alternative  to  error  diffusion. 

^This  work  was  supported  by  Office  of  Naval  Research  Grant 
N00014-93-1-0231. 


2  REVIEW  OF  DOT  DIFFUSION 

The  dot  diffusion  method  for  halftoning  has  only  one  design 
parameter,  called  class  matrix  C.  It  determines  the  order 
in  which  the  pixels  are  halftoned.  Thus,  the  pixel  positions 
(ni,n2)  of  an  image  are  divided  into  IJ  classes  according  to 
(m  mod  /,  112  mod  J)  where  I  and  J  are  constant  integers. 
Table  1  shows  an  example  of  the  class  matrix  for  /  =  J  =  8, 
used  by  Knuth.  There  are  64  class  numbers.  Let  x{ni,n2) 
be  the  contone  image  with  pixel  values  in  the  normalized 
range  [0, 1].  Starting  from  class  /c  =  1,  we  process  the  pixels 
for  increasing  values  of  k.  For  a  fixed  fc,  we  take  all  pixel 
locations  (m,  712)  belonging  to  class  k  and  define  the  halftone 
pixels  to  be 

N  if  a;(ni,7i2)  >  0.5 

x/.(ni,n2)-|Q  ifa:(ni,n2)<0.5 

We  also  define  the  error  e(ni,n2)  =  a;(ni,n2)  ”  27/1(711,712). 
We  then  look  at  the  eight  neighbors  of  (711,712)  and  replace 
the  contone  pixel  with  an  adjusted  version  for  those  neigh¬ 
bors  which  have  a  higher  class  number  (i.e.,  those  neighbors 
that  have  not  been  halftoned  yet).  To  be  specific,  neighbors 
with  higher  class  numbers  are  replaced  with 

j)  +  2e(7ii ,n2)/w  (for  orthogonal  neighbors)  (1(a)) 

x{ij)  +  e(7ii, 772)/u/  (for  diagonal  neighbors)  (1(6)) 

where  w  is  such  that  the  sum  of  errors  added  to  all  the 
neighbors  is  exactly  e(7ii,  712).  The  extra  factor  of  two  for  or¬ 
thogonal  neighbors  (i.e.,  vertically  and  horizontally  adjacent 
neighbors)  is  because  vertically  or  horizontally  oriented  error 
patterns  are  more  perceptible  than  diagonal  patterns. 

The  contone  pixels  27(711,712)  which  have  the  next  class 
number  A:  -h  1  are  then  similarly  processed.  The  pixel  values 
27(711,712)  axe  of  course  not  the  original  contone  values  but 
the  adjusted  values  according  to  earlier  diffusion  steps  (1). 
When  the  algorithm  terminates,  the  signal  (711,712)  is  the 
desired  halftone. 

This  diffusion  process  is  illustrated  in  Fig.  1.  The  num¬ 
bers  in  the  matrix  axe  elements  of  a  class  matrix  and  the 
integers  in  the  bubbles  axe  relative  weights  of  diffusion  coef¬ 
ficients.  The  neighbors  of  33  with  higher  class  number  are 
those  labeled  as  58,  45,  42,40,  63,  47.  The  error  created  at 
33  is  divided  by  the  sum  of  relative  weights  of  diffusion  co¬ 
efficients,  which  is2-i-l  +  2-|-l  +  2-{-l  =  9  in  this  case.  The 
result  of  the  division,  e,  is  the  error  to  be  diffused  to  diagonal 


Figure  1:  Error  Diffusion  from  a  point  to  the  neighbor 
points 


neighbors,  and  2e  is  diffused  to  orthogonal  neighbors.  In  64 
steps  the  algorithm  will  complete  the  halftoning  process. 

Usually  an  image  is  enhanced  [6]  before  dot  diffusion 
is  applied.  For  this  the  continuous  image  pixels  C(i,j) 
are  replaced  by  C  (i,j)  =  where  C(i,j)  = 

If.  =  0.9  the.' 

c'(i,j)  =  8C(i,j)  +  C(i,j)-  ^  C(u,v). 

0<(u-t)2  +  (t,-»2<3 

This  algorithm  is  completely  parallel  requiring  9  additions 
per  pixel,  and  no  multiplications. 

3  OPTIMIZATION  OF  CLASS  MATRIX 

Knuth  introduced  the  notion  of  barons  and  near  barons  in 
the  selection  of  his  class  matrix.  A  baron  has  only  low-class 
neighbors,  and  a  near-baron  has  one  high  class  neighbor. 
The  quantization  error  at  a  baron  cannot  be  distributed  to 
neighbors,  and  the  error  at  a  near  baron  can  be  distributed 
to  only  one  neighbor.  Knuth ’s  idea  was  that  the  number 
of  BzLTons  and  near  barons  should  therefore  be  minimized. 
He  exhibited  a  class  matrix  with  two  barons  and  two  near 
barons  (Table  1).  The  quadity  of  the  resulting  halftones  are 
still  inferior  to  error  diffusion  because  of  periodic  patterns 
similar  to  ordered  dither  methods  (see  Fig.  6).  Knuth  has 
produced  a  class  matrix  with  one  baron  and  near  baron,  but 
unfortunately  these  line  up  to  produced  objectionable  visual 
artifacts.  In  our  experience,  the  baron  near  baron  criterion 
does  not  appear  to  be  the  right  choice  for  optimization.  To 
explain  this,  define  a  A:-baron  to  be  a  position  which  has  k 
high-level  neighbors.  Thus  A:  =  0  corresponds  to  a  baron, 
fc  =  1  to  a  near  baron,  and  A:  =  8  to  an  antibaa*on.  We  have 
produced  a  class  matrix  which  minimizes  the  number  of  k- 
barons  sequentially  for  0  <  A;  <  8.  The  resulting  halftone 
quality  was  found  in  most  cases  to  be  slightly  worse  than 
Knuth ’s  original  results,  leading  us  to  conclude  that  baron 
minimization  is  not  the  right  approach.  In  Sec.  3.1  we  in¬ 
troduce  a  different  optimization  criterion  based  on  the  HVS, 
and  show  that  the  image  quality  is  significantly  improved, 
though  the  class  matrix  does  not  minimize  barons. 


3.1  Objective  Function  Based  on  Blue  Noise 
It  has  been  observed  in  the  past  that  the  error  in  a  good 
halftone  should  have  the  blue  noise  property  [3].  This 
means  that  the  noise  energy  should  mostly  be  in  the  high 
frequency  region  where  it  is  known  to  be  less  perceptible. 
We  will  show  how  to  incorporate  blue  noise  characteristics 
into  the  class  matrix  optimization. 

Imagine  that  we  have  a  constant  gray  image  x(ni,  712)  =  g 
where  0  <  5  <  1.  Let  x/i  (711,712)  denote  the  halftoned  ver¬ 
sion.  Since  the  halftone  is  supposed  to  create  the  perception 
of  the  gray  level,  the  average  number  of  dark  pixels  should  be 
equal  to  the  original  gray  level.  Typically,  therefore,  the  dark 
pixels  are  spatially  distributed  with  a  certain  average  fre¬ 
quency  fg  called  the  principal  frequency,  which  increases 
with  gray  level  g.  The  preference  for  blue  noise  [3]  (high  fre¬ 
quency  white  noise)  in  halftoning  arises  because  noise  energy 
at  a  significantly  higher  spatial  frequency  than  fg  is  not  per¬ 
ceivable.  Thus,  we  can  optimize  a  halftoning  method  for  a 
particular  gray  level  g  by  forcing  the  noise  spectrum  to  be 
concentrated  above  fg. 

This  does  not,  however,  imply  optimality  at  other  gray 
levels.  Interestingly  however,  if  the  gray  level  g  during  the 
optimization  phase  is  chosen  carefully,  the  resulting  halftones 
for  arbitrary  natural  images  are  excellent.  For  example  we 
optimized  the  class  matrix  in  the  dot  diffusion  method  for 
the  gray  level  g  —  0.0625  and  obtained  very  good  halftones 
for  natural  images  as  we  demonstrate  in  Sec.  4. 

Calculating  the  noise  spectrum.  In  order  to  implement  the 
optimization,  we  first  need  to  compute  the  noise  spectrum. 
The  halftone  pattern  x/i(7ii,  712)  for  the  gray  level  x(7ii,  712)  = 
g  has  the  error  e(7ii, 712)  =  y  —  x/i(7ii,  712),  which  is  an  iV  x  N 
image.  Imagine  that  this  is  divided  into  L  x  L  blocks  so  there 
are  B  =  {N/Vf'  blocks.  (In  our  experiment  N  =  256,  L  = 
64, H  =  16.)  Let  Em{k\.tk2)  be  the  L  x  L  DFT  of  the  the 
mth  block  of  e(7ii,  712).  We  define  the  average  noise  spectrum 
as 

1 

m=0 

Prom  this  we  compute  the  so-called  radially  averaged 
power  spectrum  Pr(A;r)  where  kr  is  a  scalar  called  the 
radial  frequency.  Since  |A:i|  and  |A;2|  range  from  0  to  L/2,  kr 
ranges  from  0  to  L/y/2.  We  take  specific  integer  values  for 
kr  and  calculate  Pr{kr)  as  follows.  For  each  chosen  kr  define 
an  annulus  A{kr)  in  the  (A;i,A;2)  plane  by  the  equation 

v/ifcf+fcf-jfcr|  <  A/2 

The  quantity  A,  which  determines  the  width  of  the  annulus 
is  chosen  as  unity  in  our  work.  With  N{kr)  denoting  the 
number  of  elements  in  A{kr),  the  radially  averaged  power 
spectrum  of  the  error  for  gray  level  g  is  then 

The  class  matrix  in  the  dot  diffusion  method  should  be  opti¬ 
mized  such  that  this  radial  spectrum  is  appropriately  shaped 
for  a  well-chosen  fixed  gray  level  g.  In  terms  of  the  radial  fre¬ 
quency  vciriable  fcr,  the  principal  frequency  for  the  halftone 
of  gray  level  g  is  given  by 

fg  =  kmaxyjg 


QUAinY  iiiCx'CCTip 


where  kmax  is  the  maximum  value  of  kr^  In  fact,  for  g  > 
0.5,  since  black  pixels  are  more  in  number,  the  halftone  is 
perceived  as  a  distribution  of  white  dots  and  we  have  to  take 
[^1  fs  ~  Umax's/!  9* 


Figure  2:  A  Typical  Figure  3:  The  Weight 
Desired  Radial  Spec-  Function  used  in  the  Op- 
trum  Characteristics  timization 

The  aim  of  the  optimization  is  to  shape  Pr(kr)  by  choice 
of  the  class  matrix  C  so  that  most  of  its  energy  is  moved 
to  the  region  kr  >  fg  (as  demonstrated  in  Fig.  2).  We 
therefore  define  the  cost  function 

rfo 

^{C,g)-  /  Pr{kr)w(kr)dkr 

Jo 

The  idea  is  to  choose  the  weighting  function  w{kr)  such  that 
upon  minimization  of  the  above  function,  Pr{kr)  has  a  low 
frequency  cutoff  at  principal  frequency  /g,  sharp  transition 
region,  and  a  flat  high  frequency  region.  The  weight  function 
was  chosen  to  be  w{kr)  =  c{kr  —  fg)^  for  0  <  kr  <  fg 
and  zero  outside.  The  constant  c  was  chosen  as  unity.  In 
the  optimization  the  integral  was  replaced  with  a  discrete 
sum.  The  choice  of  the  class  matrix  that  minimizes  this  sum 
was  performed  using  the  pairwise  exchange  algorithm 

[5]  described  below: 

1)  Randomly  order  the  numbers  in  the  class  matrix. 

2)  List  all  possible  exchanges  of  class  numbers. 

3)  If  an  exchange  does  not  reduce  cost,  restore  the  pair  to 
original  positions  and  proceed  to  the  next  pair. 

4)  If  an  exchange  does  reduce  cost,  keep  it  and  restart  the 
enumeration  from  the  beginning. 

5)  Stop  searching  if  no  further  exchanges  reduce  cost. 

6)  Repeat  the  above  steps  a  fixed  number  of  times  and  keep 
the  best  of  the  class  matrix. 

Choice  of  gray  level.  Since  the  algorithm  can  be  applied 
only  to  a  given  gray  level,  the  gray  level  should  be  chosen 
wisely,  in  order  to  get  good  halftones  for  other  gray  levels 
also.  In  our  experience  if  we  perform  this  optimization  for 
a  fixed  small  gray  level  (e.g.,  g  —  0.0625,0.1  etc.),  we  get 
good  halftones  for  natural  images  also.  Class  matrices  ob¬ 
tained  from  optimization  with  a  very  small  gray  level  will 
not  work,  because  there  is  not  much  error  to  diffuse  to  other 
points  during  the  dot  diffusion  process.  Mid  gray  levels  are 
not  suitable,  first  because  there  are  huge  diffusions  between 
points,  and  second,  even  unoptimized  algorithms  yield  per¬ 
ceptually  pleasing  halftones  for  mid  gray  anyway.  The  actual 


gray  level  used  in  the  optimization  is  ^  —  0.0625,  and  it  is 
found  experimentally.  The  class  matrix  obtained  is  shown  in 
Table  2.  Notice  that  the  optimal  class  matrix  has  several 
barons  and  near  barons. 

4  EXPERIMENTAL  RESULTS 

Fig.  4  shows  the  512  x  512  continuous  tone  peppers  im¬ 
age.  The  halftone  generated  by  using  Knuth’s  class  matrix 
is  shown  in  Fig  6,  and  the  halftone  image  generated  by  the 
optimized  class  matrix  is  shown  in  Fig  7.  It  is  clear  that  the 
new  method  is  superior  to  unoptimized  dot  diffusion  method. 
In  fact,  the  new  method  offers  a  quality  indistinguishable 
from  FS  error  diffusion  method  (Fig.  5).  Notice  that  the 
artificial  periodic  patterns  in  Fig.  6  are  absent  in  Fig.  5 
and  in  the  new  method  Fig.  7. 

5  CONCLUSION 

Even  though  dot  diffusion  offers  more  parallelism  than  error 
diffusion,  it  has  not  received  much  attention.  This  is  partly 
because  the  noise  characteristics  of  error  diffusion  method 
are  generally  regarded  as  superior.  In  this  paper  we  have 
shown  that  by  optimizing  the  class  matrix  for  blue  noise  at 
a  fixed  gray  level,  the  results  of  dot  diffusion  can  be  made 
at  least  as  pleasing  as  that  of  error  diffusion.  The  algorithm 
terminates  in  at  most  64  steps  compared  to  steps  needed 
for  error  diffusion.  Moreover,  since  the  optimized  class  ma¬ 
trix  has  eight  barons  and  eight  antibeirons  (Sec.  3.1),  the 
algorithm  can  in  fact  be  terminated  in  about  50  steps.  The 
conclusion  is  that  Knuth’s  dot  diffusion  method  with  a  care¬ 
fully  optimized  class  matrix  is  very  promising;  the  image 
quality  is  comparable  to  error  diffusion,  and  the  implemen¬ 
tation  offers  more  parallelism  than  error  diffusion. 
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Figure  4:  Original  Image,  Peppers  Figure  6:  Dot  Diffusion  with  Knuth’s  Class  Matrix 


Figure  5:  FS  Error  Diffusion 


Figure  7:  Dot  Diffusion  with  the  new  Class  Matrix 


